;+
; NAME:
;       PROCESS_GLIMPSE_IRDC
;
; PURPOSE:
;       Create the UBC_PROC postage stamp images used for the IRDC
;       distance modeling and morphological matching.
;
; CATEGORY:
;       Local file creation.
;
; CALLING SEQUENCE:
;       process_glimpse_irdc [,TONS OF DEBUGGING OPTIONS]
;
; INPUTS:
;       NONE (all inputs hardwired)
;
; OPTIONAL INPUTS:
;       GDIR      -- Directory containing the GLIMPSE Version 3 images and
;                    star-subrtracted images produced by
;                    process_glimpse_starsub.pro
;                    [Default: 'local/glimpse_v3/']. 
;       BDIR      -- Directory containing the BGPS v1.0.2 mosaics
;                    [Default: 'local/bgps/maps/']. 
;       LDIR      -- Directory containing the BGPS v1.0.2 label maps
;                    [Default: 'local/bgps/label/'].
;       SDIR      -- Directory containing the smoothed BGPS v1.0.2
;                    mosaics [Default: 'local/bgps/other/'].
;       PS_SIZE   -- Size of the postage stamp image to extract
;                    (arcmin).  [Default: 8]
;       CNUM_LIST -- List of BGPS catalog numbers to process.
;                    [Default: ALL]
;       IFACT     -- Factor multiplying MAD(irac) for histogram width
;                    about MODE(irac) to use for fitting I_MIR
;                    surface.
;
; KEYWORD PARAMETERS:
;       ADD_SUFF  -- Add the appropriate suffix to output filenames to
;                    identify the type of fit / size used.
;       FITLINEAR -- Fit a planar surface to the I_MIR rather than a
;                    quadratic.
;       VERBOSE   -- Be verbose.
;       TEST      -- Do certain testing procedures.
;
; OUTPUTS:
;       NONE (postage stamp images written to file)
;
; OPTIONAL OUTPUTS:
;       IMIR_STATS -- Statistics related to the I_MIR fitting.
;
; NOTES:
;       Requires the BGPS structure created by ASSOC_BGPS_BOLOCAT.pro
;       (IDL save file local/BGPS_v102_map_locations.sav)
;
; MODIFICATION HISTORY:
;
;       Created:  xx/xx/xx, Steve Mairs -- Initial verision, called
;                                          irdarkclouds2.pro in code
;                                          retrieved from
;                                          signals.siglab.ok.ubc.ca
;       Modified: 12/07/11, TPEB -- Modifications to allow the code to
;                                   interface with the
;                                   distance-omnibus repository, and
;                                   code cleanup for speed.
;       Modified: 12/21/11, TPEB -- Name change, and move within
;                                   repository structure.
;       Modified: 01/03/12, TPEB -- Point to new starsub files created
;                                   by process_glimpse_starsub.pro,
;                                   and other cleanup.
;       Modified: 01/11/12, TPEB -- Removed any use of the IRAC Band 1
;                                   image, which is irrelevant to this
;                                   routine.
;       Modified: 01/12/12, TPEB -- Fixed (simplified) method of
;                                   choosing which GLIMPSE or BGPS
;                                   mosaic to use, conforming with
;                                   remainder of distance-omnibus
;                                   repository.
;       Modified: 02/10/12. TPEB -- Adjusting the algorithm for
;                                   fitting the plane for I_MIR,
;                                   original code from Steve selected
;                                   pixels within 1x MAD of the mode,
;                                   but text from Erik says to use
;                                   1.47x.  Additional code changes
;                                   for testing purposes.
;       Modified: 03/26/12, TPEB -- Added new output file for the
;                                   label map.  This ensures that the
;                                   labels for morphological matching
;                                   are properly registered with the
;                                   bgps and glimpse images created
;                                   here.
;       Modified: 04/23/12, TPEB -- Added capability to fit a
;                                   quadratic surface to the I_MIR
;                                   rather than a plane.  This will be
;                                   useful for larger postage stamps.
;       Modified: 05/07/12, TPEB -- Added capacity to label filename
;                                   for postage-stamp size and type of
;                                   fit to IMIR.
;       Modified: 07/05/12, TPEB -- Modified default behavior to be
;                                   "_Q8".
;       Modified: 07/06/12, TPEB -- Updated documentation.
;       Modified: 08/07/12, TPEB -- Dealt with the 2 sources who fall
;                                   off the bottom of the BGPS mosaic.
;       Modified: 10/19/12, TPEB -- Added IRAC Band 4 scattering term
;                                   correction here so that downstream
;                                   postage stamp images may be
;                                   treated as correct.
;       Modified: 10/22/12, TPEB -- Set default mode to _q6, since the
;                                   IRAC camera FOV is 5.2' x 5.2'.
;                                   The new default size means the
;                                   scattering correction is closer to
;                                   the ideal per-IRAC-frame
;                                   correction that should have been
;                                   done within the GLIMPSE pipeline.
;                                   Also, reduced size of psf image
;                                   used to compute smoothed images.
;                                   New psf covers 10 orders of
;                                   magnitude, but no more.
;       Modified: 01/08/13, TPEB -- Changed smoothing kernel to
;                                   33" instead of 30" (BGPS group
;                                   groussing about 30").
;       Modified: 01/23/13, TPEB -- Added FITS keywords for Galactic
;                                   coordinates of the object in
;                                   question, and the estimate of the
;                                   mode from MMM.pro.
;       Modified: 01/24/13, TPEB -- Changed method for determining
;                                   "background pixels" for the I_MIR
;                                   surface fitting business.  Now
;                                   creating pixel histogram of ndata8
;                                   and selecting pixels within the
;                                   FWHM of the main peak.  In
;                                   general, this decreases the number
;                                   of pixels used in the background
;                                   fit, but the intensity range is
;                                   more closely suited to the
;                                   purpose, especially for
;                                   pathological cases.  Also added
;                                   FITS keywords related thereto.
;
;-

PRO PROCESS_GLIMPSE_IRDC, GDIR=GDIR, BDIR=BDIR, LDIR=LDIR, SDIR=sdir, $
                          PS_SIZE=ps_size, ADD_SUFF=add_suff, $
                          VERBOSE=verbose, CNUM_LIST=cnum_list, $
                          FITLINEAR=fitlinear, $
                          TEST=test, IFACT=ifact, IMIR_STATS=imir_stats
  
  ;; Parse keywords
  verbose  = KEYWORD_SET(verbose)
  silent   = ~verbose
  test     = KEYWORD_SET( test )
  fitquad  = ~ KEYWORD_SET(fitlinear)
  no_suff  = ~ KEYWORD_SET(add_suff)
  ifact = (n_elements(ifact) EQ 0) ? 1.47d : double(ifact)
  IF n_elements(ps_size) EQ 0 THEN ps_size = 6
  
  ;; Directories
  IF n_elements(GDIR) EQ 0 THEN GDIR = 'local/glimpse_v3/'
  IF n_elements(BDIR) EQ 0 THEN BDIR = 'local/bgps/maps/'
  IF n_elements(LDIR) EQ 0 THEN LDIR = 'local/bgps/label/'
  IF n_elements(SDIR) EQ 0 THEN SDIR = 'local/bgps/other/'
  
  ;; Parse optional input for size of boxes to extract around each
  ;;   source, and use default if not specified.
  ;;========
  ;; GLIMPSE: 150 pixels = 3.00'  220 pixels = 4.40'  300 pixels = 6.00'
  ;; BGPS:     16 pixels = 1.92'   25 pixels = 3.00'   33 pixels = 3.96'
  IF fitquad THEN BEGIN
     imirtype='Quadratic'
     suff = '_q'
  ENDIF ELSE BEGIN
     imirtype='Plane'
     suff = '_l'
  ENDELSE
  CASE fix(ps_size) OF
     4: BEGIN
        bdr  = 150
        bdr2 = 16
         suff += '4'
     ENDCASE
     6: BEGIN
        bdr  = 220
        bdr2 = 25
        suff += '6'
     ENDCASE
     8: BEGIN
        bdr  = 300
        bdr2 = 33
        suff += '8'
     ENDCASE
     ELSE: message,'Error: Unrecognized PS_SIZE value '+$
                   string(ps_size,format="(I0)")
  ENDCASE
  IF no_suff THEN suff = ''
  cutsize = (2.*float(bdr2)+1.)*7.2/60.
  
  ;; Begin by making a directory to place the output files
  OUTDIR = './local/glimpse_postage/'
  IF ~ FILE_TEST(OUTDIR) THEN spawn,'mkdir '+OUTDIR
  
  ;; Read in BGPS structure & parse into simple variables
  s   = read_bgps_csv(CSV=csv,/ver)
  nb  = csv.nrows
  num = s.cnum
  l   = s.glon_peak
  b   = s.glat_peak
  
  ;; Read in IDL save file with BGPS map positions
  restore,'./local/bgps_v102_map_locations.sav',/ver
  
  ;; Check for CNUM_LIST input
  IF n_elements( cnum_list ) NE 0 THEN BEGIN
     cind = WHERE_ARRAY(cnum_list, num, nb)
     IF nb EQ 0 THEN message,'Error: CNUM_LIST not matched to BGPS'
     s    = s[cind]
     num  = num[cind]
     l    = l[cind]
     b    = b[cind]
     bgps = bgps[cind]
  ENDIF
  
  ;; Next, search for the GLIMPSE/IRAC Band Images & Star-subtracted
  ;;    images, and BGPS mosaics.
  fl3 = FILE_SEARCH(GDIR+'GLM*I1.fits',              COUNT = nfl3)
  fl8 = FILE_SEARCH(GDIR+'starsub*.fits',            COUNT = nfl8)
  flb = FILE_SEARCH(BDIR+'v1.0.2_*.fits',            COUNT = nflb)
  fll = FILE_SEARCH(LDIR+'v1.0.2_*.gz',              COUNT = nfll)
  fls = FILE_SEARCH(SDIR+'*_13pca_smoothmap50.fits', COUNT = nfls)
  
  ;; Error checking
  IF nfl3 EQ 0 || nfl8 EQ 0 THEN message,'Provided GLIMPSE directory '+GDIR+$
                                         ' does not contain the required files.'
  IF nflb EQ 0 THEN message,'Provided BGPS directory '+BDIR+$
                            ' does not contain the required fields.'
  IF nflb NE nfll THEN message,'Mismatch between BGPS maps and the LABEL maps.'
  IF nfl8 NE nfl3 THEN message,'ERROR:  Mismatch between BAND 1 images and '+$
                               'STARSUB images.  May need to run '+$
                               'process_glimpse_starsub.pro.'
  
  
  
  ;; Create SMOOTHING kernel for convolving up to BGPS
  psfsize = REPLICATE( (33.0/1.2), 2)
  psf     = psf_gaussian(npixel=[155,155],fwhm=psfsize, /NORMALIZE)
  
  ;; Strings containing the name of the last file to be read in
  lastfn8 = ''
  lastfnb = ''
  
  start_t = SYSTIME(1)
  ;;===================================================================
  ;; Start looping through the objects!
  FOR p=0L, nb-1 DO BEGIN
     ;;FOR p=0L,10 DO BEGIN
     
     ellapsed = SYSTIME(1) - start_t
     frac = double(p+1)/double(nb)
     time_str = strmid(SYSTIME(0, start_t + ellapsed/frac),4,12)
     ;; Print status of progress through the catalog...
     print,string(s[p].cnum,nb,float(p+1)/float(nb)*100.,time_str,format=$
                  "('Running Bolocat object #',I4,' / ',I4,'   (',F5.1,'%)"+$
                  "  ',A0)")
     
     time0 = systime(1)
     
     ;; Choose which GLIMPSE file to use
     glon = string( round(s[p].glon_peak / 3.) * 3., format="(I03)")
     IF glon EQ '360' THEN glon = '000'
     thisfn8 = GDIR + 'starsub.'+glon+'00.I4.fits'
     
     ;; Choose which BGPS file to use
     thisfnb = BDIR + bgps[p].filename
     
     index  = WHERE(fl8 EQ thisfn8, nind1)
     index2 = WHERE(flb EQ thisfnb, nind2)
     
     ;; If out of range of GLIMPSE, skip this source
     IF verbose THEN print,'NINDs for IRAC & BGPS: ',nind1,nind2
     IF nind1 EQ 0 || nind2 EQ 0 THEN CONTINUE
     
     ;;=================================================================
     ;; Optimization - trying to get the script to run faster, i.e. we
     ;;   do not need to read in a new fits file each time if it is
     ;;   indeed the same file
     
     IF thisfn8 NE lastfn8 THEN BEGIN
        message,'Reading in new GLIMPSE mosaic... '+thisfn8,/inf
        lastfn8 = thisfn8
        undefine,data8,adata8,hd8,ahd8,astrom8
        
        ;; Read in file & perform header astrometery (force FLOAT)
        data8 = TEMPORARY(float(readfits(thisfn8,hd8)))
        extast,hd8,astrom8
        new8 = 1b
        
        ;; Muck with the FITS header
        sxaddpar,hd8,'BITPIX',-32
        sxaddhist, '==============================', hd8, /COMMENT
        hist = ['PROCESS_GLIMPSE_IRDC: '+systime(0),$
                'Created IRDC postage stamp images']
        sxaddhist, hist, hd8
        sxaddpar, hd8, 'SSUBFILE', sxpar(hd8,'FILENAME'), $
                  ' Star-subtracted GLIMPSE image'
        sxaddpar, hd8, 'CUTSIZE', cutsize, ' [arcmin] Size of postage stamp', $
                  FORMAT="F0.2"
        sxdelpar, hd8, ['DELTA-X','DELTA-Y','BORDER','RA','DEC']
        FOR rmi = 1, 300 DO BEGIN
           sxdelpar, hd8, string(rmi,format="('AOR',I03)")
           sxdelpar, hd8, string(rmi,format="('DSID',I03)")
        ENDFOR
     ENDIF      
     
     ;; Now, do BGPS image
     IF thisfnb NE lastfnb THEN BEGIN
        message,'Reading in new BGPS mosaic... '+thisfnb,/inf
        lastfnb = thisfnb
        undefine,datab,hdb,astromb
        
        ;; Read in file & perform header astrometery  (force FLOAT)
        datab = TEMPORARY(float(readfits(thisfnb,hdb)))
        label = TEMPORARY(fix(readfits(LDIR + bgps[p].labelname,hdl)))
        datas = TEMPORARY(float(readfits(SDIR + bgps[p].smooname,hds)))
        extast,hdb,astromb
        
        ;; Muck with the FITS header
        basename = strmid(thisfnb,strpos(thisfnb,'/',/reverse_search)+1)
        sxaddpar,hdb,'BITPIX',-32
        sxaddpar,hdb,'FILENAME',basename,' Name of this file',$
                 AFTER='WAVELENG'
        sxaddhist, '==============================', hdb, /COMMENT
        hist = ['PROCESS_GLIMPSE_IRDC: '+systime(0),$
                'Created IRDC postage stamp images']
        sxaddhist, hist, hdb
        sxaddpar,hdb,'ORIGFILE',sxpar(hdb,'FILENAME'),' Original BGPS mosaic'
        sxaddpar, hdb, 'CUTSIZE', cutsize, ' [arcmin] Size of postage stamp', $
                  FORMAT="F0.2"
        FOR rmi = 1, 300 DO BEGIN
           sxdelpar, hdb, string(rmi,format="('AOR',I03)")
           sxdelpar, hdb, string(rmi,format="('DSID',I03)")
        ENDFOR
     ENDIF 
     
     ;; Perform header astrometry
     ad2xy,l[p],b[p],astrom8,x8,y8
     ad2xy,l[p],b[p],astromb,xb,yb
     
     ;; Define the boxes for postage stamp image extraction, rounding
     ;;   to nearest pixel
     sz8 = [astrom8.naxis[0],astrom8.naxis[1]]
     szb = [astromb.naxis[0],astromb.naxis[1]]
     
     box     = long(round([(x8-bdr)>0,(x8+bdr)<(sz8[0]-1),$
                           (y8-bdr)>0,(y8+bdr)<(sz8[1]-1)]))
     bgpsbox = long(round([(xb-bdr2)>0,(xb+bdr2)<(szb[0]-1),$
                           (yb-bdr2)>0,(yb+bdr2)<(szb[1]-1)]))
     
     ;; Go about extracting the proper sections of the images.
     hextract,data8,hd8,ndata8,nhd8,box[0],box[1],box[2],box[3],SILENT=silent
     hextract,datab,hdb,bgpsstamp,bgpsstamphd,$
              bgpsbox[0],bgpsbox[1],bgpsbox[2],bgpsbox[3],SILENT=silent
     hextract,label,hdl,lblstamp,lblstamphd,$
              bgpsbox[0],bgpsbox[1],bgpsbox[2],bgpsbox[3],SILENT=silent
     hextract,datas,hds,smoostamp,smoostamphd,$
              bgpsbox[0],bgpsbox[1],bgpsbox[2],bgpsbox[3],SILENT=silent
     
     ;; Add GLON / GLAT for this source to the FITS header
     sxaddpar,nhd8,'GLON',l[p],' [deg] Galactic Longitude for this object',$
              FORMAT="F0.4"
     sxaddpar,nhd8,'GLAT',b[p],' [deg] Galactic Latitude for this object',$
              FORMAT="F0.4"
     
     ;; Check on the 2 sources falling off the bottom of a BGPS mosaic
     IF (yb-bdr2) LT 0 THEN BEGIN
        below = long(round(abs(yb-bdr2)))
        newstmp = fltarr(2*bdr2+1,2*bdr2+1)
        newstmp[*,below:*] = bgpsstamp
        sxaddpar,bgpsstamphd,'NAXIS2',2*bdr2+1,format="(I0)"
        sxaddpar,bgpsstamphd,'CRPIX2',sxpar(bgpsstamphd,'CRPIX2')+below
        bgpsstamp = TEMPORARY(newstmp)
     ENDIF
     
     ;;================================================================
     ;; Section of code dealing with fitting the I_MIR to the cutout
     
     boxsz = size(ndata8,/DIM)
     xarr = cmreplicate(findgen(boxsz[0]),boxsz[1]) ; The column array
     yarr = transpose(cmreplicate(findgen(boxsz[1]),boxsz[0])) ; The row array
     
     ;; Estimate the sky background in a stellar contaminated field.
     mmm,ndata8,mo,skysig,/SILENT
     sxaddpar,nhd8,'SKYMODE',mo,' [MJy/sr] Mode of image from MMM.pro',$
              FORMAT="F0.3"
     sxaddpar,nhd8,'SKYSIG',skysig,FORMAT="F0.3",$
              ' [MJy/sr] Sky sigma from MMM.pro (-1 = FAIL)'
     IF skysig EQ -1 THEN BEGIN
        ;; I guess, at one time, I thought the program should do
        ;;   something here?
        
     ENDIF
     xcen = (box[1]-box[0])/2
     ycen = (box[3]-box[2])/2
     
     ;; Maybe need to divide second term in where function by 2? --
     ;;   IFACT = 1.47 from paragraph Erik sent for IRDC paper
     index = where(abs(ndata8-mo) LT mad(ndata8)*ifact)
     
     ;;============
     ;; Alternative version: Make pixel histogram of HASTROM'd
     ;;   image with 1 MJy/sr bins, and take all pixels from
     ;;   the 1.2"-resolution image w/in FWHM
     histbin = 1.0
     hastrom,ndata8,nhd8,test8,thd8,bgpsstamphd,missing=0.
     plothist,/noplot,test8,Iarr,Narr,bin=histbin
     ihist = where(Narr GE 0.5*max(Narr))
     fitr = minmax(Iarr[ihist]) + ([-1.,1.]*histbin/2.)
     index = where(ndata8 GE fitr[0] AND ndata8 LE fitr[1], nindex)
     sxaddpar,nhd8,'IFIT_MIN',fitr[0],FORMAT="F0.3",$
              ' [MJy/sr] Minimum pix value used for I_MIR fit'
     sxaddpar,nhd8,'IFIT_MAX',fitr[1],FORMAT="F0.3",$
              ' [MJy/sr] Maximum pix value used for I_MIR fit'
     sxaddpar,nhd8,'IFIT_NPIX',nindex,' Number of pixels used for the I_MIR fit'
     ;;============
     
     IF test THEN BEGIN
        mask        = ndata8*0.
        mask[index] = 1.
        plotmap,mask,nhd8,title='IFACT = '+string(ifact,format="(F0.1)")
     ENDIF
     
     ;;=================================================================
     ;; Perform the **LEAST SQUARES FIT**
     ;; mat ## fitvec = RHS
     
     ;; (x,y,z) coordinates of the pixels used for the fit     
     xp = xarr[index]-xcen
     yp = yarr[index]-ycen
     zp = ndata8[index]
     
     ;; Calculate the TERMS used for the matrix
     nterm = fitquad ? 6 : 3
     mat = dblarr(nterm,nterm)
     
     mat[0,0] = n_elements(xp)
     mat[1,0] = total(xp)
     mat[2,0] = total(yp)
     mat[1,1] = total(xp*xp)
     mat[2,1] = total(xp*yp)
     mat[2,2] = total(yp*yp)
     IF fitquad THEN BEGIN
        mat[3,0] = total(xp*xp)
        mat[4,0] = total(yp*yp)
        mat[5,0] = total(xp*yp)
        mat[3,1] = total(xp*xp*xp)
        mat[4,1] = total(xp*yp*yp)
        mat[5,1] = total(xp*xp*yp)
        mat[3,2] = total(xp*xp*yp)
        mat[4,2] = total(yp*yp*yp)
        mat[5,2] = total(xp*yp*yp)
        mat[3,3] = total(xp*xp*xp*xp)
        mat[4,3] = total(xp*xp*yp*yp)
        mat[5,3] = total(xp*xp*xp*yp)
        mat[4,4] = total(yp*yp*yp*yp)
        mat[5,4] = total(xp*yp*yp*yp)
        mat[5,5] = total(xp*xp*yp*yp)
     ENDIF
     
     ;; Fill in the symmetric matrix
     FOR ii=0,nterm-1 DO mat[ii,*] = mat[*,ii]
     
     ;; The right hand side of the matrix equation
     RHS    = dblarr(1,nterm)
     RHS[0] = total(zp)
     RHS[1] = total(xp*zp)
     RHS[2] = total(yp*zp)
     IF fitquad THEN BEGIN
        RHS[3] = total(xp*xp*zp)
        RHS[4] = total(yp*yp*zp)
        RHS[5] = total(xp*yp*zp)
     ENDIF    
     
     ;; Calculate the fit vector
     fitvec = (invert(mat))##RHS
     
     ;; Calculate I_MIR
     Imir = fitvec[0] + fitvec[1]*(xarr-xcen) + fitvec[2]*(yarr-ycen)
     IF fitquad THEN $
        Imir += fitvec[3]*(xarr-xcen)*(xarr-xcen) + $
                fitvec[4]*(yarr-ycen)*(yarr-ycen) + $
                fitvec[5]*(xarr-xcen)*(yarr-ycen)
     
     IF test THEN $
        imir_stats = m4_stat(Imir)
     
     ;; Clear memory -- Cleaner in the code this way than tons of
     ;;                 TEMPORARY() statements.
     undefine,index,xp,yp,zp,xarr,yarr
     
     
     ;;============================================================
     ;; IRAC Band 4 suffers from pretty bad scattered light within the
     ;; camera.  No light is lost, but it is internally reflected all
     ;; over the place.  Therefore, there is a photometric aperture
     ;; correction that must be applied to surface brightness
     ;; measurements.  For extended diffuse emission (read: Galactic
     ;; Plane), this correction factor is 0.737, meaning the true
     ;; surface brightness is only 73.7% of that in the GLIMPSE
     ;; mosaics.  Furthermore, for IRDCs the dark region is filled
     ;; with scattered BRIGHT emission.  The propagation of this
     ;; correction leads to the TRUE CONTRAST = C_meas(1 + \xi), where
     ;; \xi is 1 - correction, or 0.263.
     xi = 0.263
     
     ;; We deal with the scattering here by subtracting X = xi * IMIR
     ;; from both the IMIR fit (above) and the ndata8 cutout.  Add
     ;; value of X to the FITS headers.
     
     ;; Xscat is xi times the mean of IMIR rather than a
     ;; position-dependent array:
     Xscat = xi * mean(Imir)
     
     Imir -= Xscat
     ndata8 -= Xscat
     
     sxaddpar,nhd8,'X_SCAT',Xscat,$
              ' [MJy/sr] Scattered intensity subtracted',FORMAT="F0.3"
     
     
     ;;================================================================
     ;; Section of code dealing with smoothing
     
     ;; Median smooth regular Band 4 image to remove effects of
     ;;    badly-subtracted stars
     ndata8 = median(ndata8 > 0., 5)
     
     ;; Make smoothed image
     smooth30 = float(convolve(ndata8,psf,FT_PSF=psf_ft,/NO_PAD) > 0.)
     
     print,'After convolve:  ',systime(1) - time0
     
     ;;================================================================
     ;; Start the great FITS-writing!
     ;; NOTE: resampling up to BGPS pixel scale done at this point!!!
     
     ;; Non-smoothed image
     hastrom,ndata8,nhd8,ndata8,svhd8,bgpsstamphd,missing=0.
     sxaddpar,svhd8,'PIXSCAL1',abs(sxpar(svhd8,'CD1_1'))*3600.,format="F0.3"
     sxaddpar,svhd8,'PIXSCAL2',abs(sxpar(svhd8,'CD2_2'))*3600.,format="F0.3"
     wfn = 'nonsmooth'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,svhd8,'FILENAME',wfn
     writefits,OUTDIR+wfn,TEMPORARY(float(TEMPORARY(ndata8))),TEMPORARY(svhd8)
     
     ;; Smoothed image
     hastrom,smooth30,nhd8,smooth30,svhd8,bgpsstamphd,missing=0.
     sxaddpar,svhd8,'PIXSCAL1',abs(sxpar(svhd8,'CD1_1'))*3600.,format="F0.3"
     sxaddpar,svhd8,'PIXSCAL2',abs(sxpar(svhd8,'CD2_2'))*3600.,format="F0.3"
     wfn = 'smooth33arc'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,svhd8,'FILENAME',wfn
     writefits,OUTDIR+wfn,TEMPORARY(float(TEMPORARY(smooth30))),TEMPORARY(svhd8)
     
     ;; Imir image
     hastrom,Imir,nhd8,Imir,svhd8,bgpsstamphd,missing=0.
     sxaddpar,svhd8,'PIXSCAL1',abs(sxpar(svhd8,'CD1_1'))*3600.,format="F0.3"
     sxaddpar,svhd8,'PIXSCAL2',abs(sxpar(svhd8,'CD2_2'))*3600.,format="F0.3"
     wfn = 'Imir'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,svhd8,'FILENAME',wfn
     sxaddpar,svhd8,'IFACT',ifact,' Times MAD around mode to fit IMIR',$
              FORMAT="F0.2"
     sxaddpar,svhd8,'IMIRTYPE',imirtype,' Type of fit to I_MIR'
     writefits,OUTDIR+wfn,TEMPORARY(float(TEMPORARY(Imir))),TEMPORARY(svhd8)
     undefine,nhd8
     
     ;; BGPS label
     hastrom,lblstamp,lblstamphd,bgpsstamphd,missing=0.
     wfn = 'label'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,lblstamphd,'FILENAME',wfn
     sxaddpar,lblstamphd,'BITPIX',8
     writefits,OUTDIR+wfn,TEMPORARY(byte(TEMPORARY(lblstamp) EQ num[p])),$
               TEMPORARY(lblstamphd)
     
     ;; BGPS smoothed
     hastrom,smoostamp,smoostamphd,bgpsstamphd,missing=0.
     wfn = 'bsmoo'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,smoostamphd,'FILENAME',wfn
     writefits,OUTDIR+wfn,TEMPORARY(float(TEMPORARY(smoostamp))),$
               TEMPORARY(smoostamphd)
     
     ;; BGPS image
     wfn = 'bgps'+string(num[p],format="(I04)")+suff+'.fits'
     sxaddpar,bgpsstamphd,'FILENAME',wfn
     writefits,OUTDIR+wfn,TEMPORARY(float(TEMPORARY(bgpsstamp))),$
               TEMPORARY(bgpsstamphd)
     
  ENDFOR   ;; END of source looping
  
  ;; Clear out memory
  undefine,data8,hd8,astrom8
  undefine,datab,hdb,astromb,label,hdl,datas,hds
  
  ;;======================================
  ;; End of Code
  print,'End Time:' & spawn,'date'
  
END
